The present application hereby claims priority under 35 U.S.C. xc2xa7119 on German patent publication number DE 10155089.8 filed Nov. 9, 2001, the entire contents of which are hereby incorporated herein by reference.
Image data of an examined measurement object can be obtained using modern medical diagnosis methods, such as X-ray computed tomography (CT), for example. The examined measurement object is generally a patient.
X-ray computed tomographyxe2x80x94abbreviated to CT hereinafterxe2x80x94is a special X-ray recording method which fundamentally differs from the traditional X-ray tomography method in terms of image construction. CT recordings yield transverse sectional images, that is to say images of body layers which are oriented essentially perpendicularly to the body axis. The tissue-specific physical quantity represented in the image is the distribution of the attenuation value of X-ray radiation xcexc(x, y) in the sectional plane. The CT image is obtained by reconstructing the one-dimensional projections of the two-dimensional distribution of xcexc(x, y), which projections are supplied by the measurement system used, from numerous different viewing angles.
The projection data are determined from the intensity I of an X-ray after it has passed through the layer to be imaged and its original intensity I0 at the X-ray source in accordance with the absorption law                               ln          ⁢                      I                          I              0                                      =                  -                                    ∫              L                        ⁢                                          μ                ⁡                                  (                                      x                    ,                    y                                    )                                            ⁢                              ⅆ                l                                                                        (        1        )            
The integration distance L represents the path of the X-ray considered through the two-dimensional attenuation distribution xcexc(x, y). An image projection is then composed of the measured values of the line integrals through the object layer, with the measured values being obtained with the X-rays of a viewing direction.
The projectionsxe2x80x94characterized by the projection angle xcex1xe2x80x94originating from many different directions are obtained by a combined X-ray tube detector system which rotates about the object in the layer plane. The most common apparatuses at the present time are so-called xe2x80x9cfan beam apparatusesxe2x80x9d, in which a tube and an array of detectors (a linear arrangement of detectors) rotate in the layer plane jointly about a center of rotation, which is also the center of the circular measurement field. xe2x80x9cParallel beam apparatusesxe2x80x9d, which exhibit very long measurement times, are not explained here. It shall be pointed out, however, that a transformation from fan to parallel projections and vice versa is possible, so that the present invention, which will be explained with reference to a fan beam apparatus, can also be used for parallel beam apparatuses without any restriction.
In the case of fan beam geometry, a CT recording comprises line integral measured valuesxe2x80x94ln(I/I0) of arriving beams which are characterized by a two-dimensional linkage of the projection angle xcex1xcex5[0,2xcfx80) and the fan angles xcex2xcex5[xe2x88x92xcex20,xcex20] defining the detector positions (xcex20 is half the fan aperture angle). Since the measurement system only has a finite number k of detector elements and a measurement comprises a finite number y of projections, this linkage is discrete and can be represented by a matrix:
{tilde over (p)}(xcex1y,xcex2k)[0,2xcfx80)xc3x97[xe2x88x92xcex20,xcex20]xe2x80x83xe2x80x83(2)
or 
{tilde over (p)}(y,k): (1,2, . . . NP)xc3x97(1,2, . . . NS)xe2x80x83xe2x80x83(3)
The matrix {tilde over (p)}(y, k) is called a sinogram for fan beam geometry. The projection number y and the channel number k are of the order of magnitude of 1000.
The principle of image reconstruction in CT by calculating the xcexc value distribution will not be discussed further, for the sake of brevity. This is illustrated in detail, for example, in xe2x80x9cBildgebende Systeme fxc3xcr die medizinische Diagnostikxe2x80x9d [xe2x80x9cImaging Systems for Medical Diagnosisxe2x80x9d], third edition, Munich: Publicis MCD Verlag, 1995, ed.: Morneburg, Heinz, ISBN 3-89578-002-2.
However, the task of image reconstruction is not yet concluded with the calculation of the xcexc value distribution of the radiographed layer. In the medical field of application, the distribution of the attenuation coefficient xcexc represents only an anatomical structure which still has to be represented in the form of an X-ray image.
After a proposal by G. N. Hounsfield, it has become generally customary to transform the values of the linear attenuation coefficient xcexc (which has the dimension unit cmxe2x88x921) to a dimensionless scale in which water acquires the value 0 and air the value xe2x88x921000. The conversion formula to this xe2x80x9cCT numberxe2x80x9d reads as follows:                               CT          ⁢                      xe2x80x83                    ⁢          number                =                                            μ              -                              μ                Water                                                    μ              Water                                ⁢          1000                                    (        4        )            
The unit of the CT number is called a xe2x80x9cHounsfield Unitxe2x80x9d (HU). This scale is highly suited to the representation of anatomical tissue since the unit HU expresses the deviation in thousandths of xcexcwater and the xcexc values of most substances inherent to the body differ only little from the xcexc value of water. From the range of numbers (from xe2x88x921000 for air to approximately 3000), only integers are used as carriers of image information.
However, the representation of the entire scale range of about 4000 values would far surpass the discrimination capability of the human eye. Moreover, the observer is often interested only in a small excerpt from the range of attenuation values, e.g. the differentiation of grey and white brain substance, which differ only by about 10 HU.
For this reason, so-called image windowing is used. In this case, only part of the CT value scale is selected and spread over all available grey shades. Even small attenuation differences within the chosen window thus become perceptible grey tone differences, while all the CT values below the window are represented black and all the CT values above the window white. The image window can be varied arbitrarily but in terms of its central level and in terms of its width.
The image data obtained usually contained not only the desired image information of the examined measurement object but also information which can be attributed to disturbing influences during the measurement operation.
Generally, a distinction is made between two different categories of problems which reduce the quality of the image data obtained: image noise and artifacts. These two problems will be explained in more detail below.
Image noise can in turn be subdivided into a plurality of causes.
The main part of the image noise is brought about by the quantum noise which results from the fact that each radiation comprises a finite number of quanta, so that the number of measured quanta always virtuates normally distributed about an average value.
Further causes of the image noise are the usually not exactly monochromatic quanta of the X-ray tubes that can be realized in practice, and scattered radiation due to interactions between the X-ray radiation used and the electron shell of atoms during transmission through the measurement object.
Artifacts are also subdivided further: aliasing, partial volume artifacts, hardening artifacts and motion artifacts are typical artifacts whose occurrence depends in particular on the geometry or a motion of the measurement object.
Effects corresponding to the above-described image noise and the artifacts can also be found in other imaging systems for medical diagnosis.
Ring artifacts constitute a particular form of artifacts the cause of which is to be sought primarily in the computer tomography imaging system itself that is used:
As already described above, a plurality of detectors (up to 1000) are used in fan beam apparatuses. Therefore, there is the possibility of inadequate calibration of the individual detectors. In other words, identical attenuations of the radiation penetrating through the measurement object are measured differently by different detectors.
In the case of inadequate calibration of the individual detectors of a computer tomograph, the image data obtained have concentric rings or partial ring arcs about the center of rotation, which have no actual reference to the measurement object considered. This occurs on account of the rotation of the beam source and of the detectors about the measurement object during the measurement operation. Such disturbances in the image data are called ring artifacts.
The order of magnitude of the differential channel errors is approximately a factor of less than or equal to xc2x12xc3x9710xe2x88x923 of the detected intensity. Although these errors cause the measured value to deviate from the xe2x80x9ctruexe2x80x9d value usually only by a few thousandths, they bring about distinctly visible ring artifacts in the image. The rings or partial rings have an amplitude of xc2x120 HU.
In order to eliminate ring artifacts in the sinogram or in the X-ray image, there are various procedures in the prior art:
The patent specification U.S. Pat. No. 5,210,688 A illustrates a method for suppressing apparatus-dictated xe2x80x9cdiscontinuitiesxe2x80x9d in the sinogram which are presented as ring artifacts in the later CT image. The method is based on the determination and subsequent subtraction of a correction sinogram from the starting sinogram. To that end, the data are fourier-transformed in order then to separate channel errors from object structures in the frequency domain by low-pass filtering.
In the patent specification U.S. Pat. No. 5,745,542 A, actual channel errors are separated from signal structures in the correction sinogram by evaluation of a histogram, low-pass filtering being effected in the projection direction and high-pass filtering being effected in the channel direction.
In the patent specification U.S. Pat. No. 6,115,445 A, in order to create the correction sinogram, likewise a low-pass filtering is carried out in the projection direction and a high-pass filtering is carried out in the channel direction. However, the differentiation between useable object signals and error signals to be eliminated is effected by means of a weighting and limiting unit which carries out an error signal estimation using so-called weighting factors.
The patent specification U.S. Pat. No. 6,094,467 A discloses, for the purpose of identifying object edges (in the case of high-intensity metal implants with respect to the CT imaging, which as such considerably impair the CT image quality), evaluating A(i,xcex8) through the second derivative of the attenuation values (sinogram values) according to the channel number i (xcex8 is the respective projection).
Further methods for eliminating ring artifacts in the sinogram or in the X-ray image which comprise parts of the abovementioned method steps are the algorithms known:
a) Raw data balancing method
In this case, the channel errors which are identified in the sinogram by line-like structures are detected and corrected.
b) Image balancing method
Rings and partial rings represented in the image on account of channel errors are identified and corrected directly in the image.
It has been shown empirically that b) generally acts more efficiently than a), but in the vicinity of the center of rotation the correction can itself produce artifacts on account of a lack of statistics during the detection of rings. Therefore, the raw data balancing method is advantageous for the channels in the surroundings of the center of rotationxe2x80x94the term channel is used to denote the connecting straight line from the (point) X-ray source to a detector element.
The present application, in a preferred embodiment, can constitute a considerable further development of the raw data balancing method, for which reason the latter will be described below:
The method may be based on the following assumptions:
1. The ring artifacts to be corrected ought to sweep over a minimum angle of 30xc2x0 in the image. Accordingly, a channel error lasts for at least Np/12 successive projections.
2. The error amplitude (the ring) must not exceed a certain limit value (presently 15 attenuation units).
The so-called xe2x80x9cbalancing procedurexe2x80x9d servesxe2x80x94as mentioned abovexe2x80x94for correcting CT raw data (which are represented in the sinogram) whose associated image has ring-shaped artifacts. These errors must be identified by the method and distinguished from xe2x80x9cgenuinexe2x80x9d structures of similar appearance.
The balancing method has the following substeps:
I) Compression
xe2x80x83Firstly, by use of an averaging compression in xcex1, the volume of data is reduced and noise smoothing is carried out.
II) High-pass filtering
xe2x80x83A high-pass filtering in the direction xcex2 subsequently emphasises the differential channel errors.
III) xcex1 low-pass filter 1
xe2x80x83An xcex1 low-pass filter attenuates short-range structures, which makes it possible to distinguish channel errors from high-frequency structures.
IV) xcex1 differentiator
xe2x80x83Through xcex1 difference formation, signal structures not oriented in the xcex1 direction remain discernible.
V) Decision unit
xe2x80x83A decision as to whether or not there are channel errors is taken with the aid of an amplitude and an xcex1 gradient threshold value decision.
VI) xcex1 low-pass filter 2
xe2x80x83The short-range structures in xcex1 which are present after the decision and are not channel errors are weakened in terms of their influence by a second xcex1 low-pass filter.
VII) Decompression
VIII) Correction
xe2x80x83The correction is finally effected by subtracting the correction sinogram obtained from the measured sinogram.
The respective steps will now be specifically discussed in more detail.
In this case, reference is made to the representation of the sinogram in which the preprocessed CT valuesxe2x80x94called CT raw dataxe2x80x94represent the attenuation values S(y, k). The channels k of a projection are plotted horizontally and the projection numbers yxe2x80x94both beginning with 1 in each casexe2x80x94are plotted vertically:
I) Compression
In order to reduce the computation time and because the channel errors to be removed have a minimum extent in the projection direction, it is expedient, for the further operations, to reduce the number of projections used for the correction. This is done by arithmetic averaging of N projections, i.e. a new sinogram is created which now contains NPRO/N projections:                               Comp          ⁡                      (                          y              ,              k                        )                          =                              (                          1              N                        )                    ⁢                                    ∑                              i                =                                                      N                    ·                                          (                                              y                        -                        1                                            )                                                        +                  1                                                            y                ·                N                                      ⁢                          S              ⁡                              (                                  i                  ,                  k                                )                                                                        (        5        )            
where y=1,2 . . . NPRO/N and NPRO/N=NCOMP.
The factor N depends on the number of projections per circulation. N is chosen such that the distance between two compressed projections corresponds to a circulation angle xcex1 of approximately 3xc2x0. These are approximately 120 projections in the case of a full circulation.
As stated, although a compression is expedient, it does not necessarily have to be carried out.
II) High-pass Filtering
In compressed sinograms, too, the channel errors continue in successive compressed projections. Since the errors are very small, however, and, moreover, only the differential amplitude error with regard to the adjacent channels is intended to be corrected, it is possible to perform a suitable high-pass filtering in the channel direction (xcex2 direction):                               High          ⁢                      -                    ⁢                      pass            ⁡                          (                              y                ,                k                            )                                      =                              ∑                          i              =                              L                2                                                    L              2                                ⁢                                    Comp              ⁡                              (                                  y                  ,                                      k                    -                    i                                                  )                                      ·                          hp              ⁡                              (                i                )                                                                        (        6        )            
The effect of the high-pass filter in the image and sinogram region is characterized in that low-frequency components are eliminated but ring artifacts and structures with high-frequency components in the radial direction are emphasized. Since the object has a long-wave profile, it is filtered out by the high-pass filter. What remain, on the one hand, are high-frequency channel errors which are removed at the end of the method by subtraction of the correction sinogram from the original sinogram. On the other hand, sharp object edges are also identified as high-frequency structures, not filtered out and included in the correction sinogram. Thus, important object structures would be incorrectly eliminated after the sinogram correction. This is to be prevented, however, in the xcex1 differentiator in step IV).
In order to eliminate so-called outliers, a median filter may likewise optionally be used after the high-pass filtering. The median filter also avoids the situation where channel errors are easily smeared out in the high-pass filtering.
III) xcex1 Low-pass Filter 1
As mentioned in the introduction, ring artifacts in the image ought to have a minimum extent of 30 degrees. xe2x80x9cGenuinexe2x80x9d structures which are arranged in a ring-shaped manner but have a smaller extent are not to be removed. To that end, a smoothing operation is introduced in the projection direction, which smoothing operation reduces the amplitude of short structures but retains the amplitude of extended rings.                               Low          ⁢                      -                    ⁢                      pass1            ⁡                          (                              y                ,                k                            )                                      =                              ∑                          i              =                              -                                  L                  2                                                                    L              2                                ⁢                      High            ⁢                          -                        ⁢                                          pass                ⁡                                  (                                                            y                      -                      i                                        ,                    k                                    )                                            ·                              tp                ⁡                                  (                  i                  )                                                                                        (        7        )            
In addition, this reduces the noise and thus increases the deduction reliability. This smoothing in the projection direction makes it possible to isolate channel errors and eliminate noise.
It may be, however, that objects (e.g. skull bones) cause lines in the sinogram which are very similar to those of the channel or amplitude errors. Such structures must likewise be eliminated in order that they are not retouched away during the correction. To that end, the step of differentiation in the projection direction is used.
IV) xcex1 Differentiator
In order to detect structures that are not (exactly) ring-shaped and are distinguished by great change in the projection direction, the smoothed signal is xe2x80x9cdifferentiatedxe2x80x9d in the projection direction.
Diff(y,k)=Low-pass1(y,k)xe2x88x92Low-pass1(yxe2x88x921,k)xe2x80x83xe2x80x83(8)
This operation is an approximation of the first derivatives in the xcex1 direction or of the subtractive parts of the first derivative (gradient) in the image region:                                           ∂                          ∂              α                                ⁢          Low          ⁢                      -                    ⁢                      pass1            ⁡                          (                              y                ,                k                            )                                      =                                            Low              ⁢                              -                            ⁢                              pass1                ⁡                                  (                                                            y                      +                                              Δ                        ⁢                                                  xe2x80x83                                                ⁢                        y                                                              ,                    k                                    )                                                      -                          Low              ⁢                              -                            ⁢                              pass1                ⁡                                  (                                      y                    ,                    k                                    )                                                                          Δ            ⁢                          xe2x80x83                        ⁢            y                                              (        9        )            
If in equation (9) xcex94y=1 (best possible approximation), the following results                                                         ∂                              xe2x80x83                                                    ∂              α                                ⁢          Low          ⁢                      -                    ⁢                      pass1            ⁡                          (                              y                ,                k                            )                                      ≈                  Diff          ⁡                      (                          y              ,              k                        )                                              (        10        )            
Although signal structures that are not oriented in the xcex1 direction and signal structures that are oriented in the xcex1 direction with rapidly variable amplitude are altered, they nonetheless remain discernible. Data structures of whole rings and partial rings are eliminated. In other words, structures which cannot be assigned to a channel error, but rather are caused by absorbing objects, are eliminated, to be precise on the basis of the assumption that a channel error varies only slightly over time, but an object edge effects a rapid variation.
V) Decision Unit
The decision as to whether or not rings are present is made through two amplitude thresholds. Specifically, it may be that a structure has been detected which, incorrectly, was not a ring but rather, however for example, a sharp bone edge (such a structure generally has a very high amplitude). If the magnitude of the derivative exceeds one of the thresholds defined below, then the assigned data point is removed in order that this in actual fact anatomical structure is not taken into account later during the sinogram subtraction.
1) If the magnitude of the absolute value of the differentiated signal Diff(y, k) lies below the gradient threshold S4 less than 0, it is assumed that what is involved is a ring whose amplitude has been calculated in the sinogram Low-pass1(y, k). Otherwise, this value is set to 0.
2) If the magnitude of the absolute value of the sinogram |Low-pass1(y,k)| is greater than an amplitude threshold S3 greater than 0, then the value is limited to xc2x1S3. This is intended to reduce the effect of incorrect corrections.
The threshold value operations can be described as follows:
The threshold S3 is channel-independent, and S4 is channel-dependent. The threshold S4i holds true in a central channel region, and S4a holds true peripherally (S4i greater than S4 a). Identification conditions in the outer channel region are stricter than in the central channel region (S4i=2*S4a holds true in the outer region). The inner region of a projection usually contains higher-frequency data structures than the outer region. The risk of elimination of channel errors through an excessively low gradient threshold is then greater.
To summarize, it may be stated that, assuming that channel errors have amplitudes in a restricted range and, as a result, the data are limited by a threshold, a generation of artifacts on account of objects not identified in step IV) can be counteracted by step V).
VI) xcex1 Low-pass Filter 2
A concluding xcex1 low-pass filtering of the sinogram generated by the decision unit is intended to smooth the sudden amplitude changes generated by the threshold value operations.
The edge effects discussed above are thus largely eliminated. Since a minimum extent of 30xc2x0 was assumed in the ring detection, object structures having an extent of a few millimeters would also be identified as rings and removed in the detector center (in the central channel region). Therefore, a smoothing of significantly greater range by means of a second xcex1 low-pass filter is used in a region around the detector center. For the central region of the detector channels k, equation (12) holds true where e.g. L=25:                              Low          ⁢                      -                    ⁢                      pass2            ⁡                          (                              y                ,                k                            )                                      =                              ∑                          i              =                              -                                  L                  2                                                                    L              2                                ⁢                                    Dec              ⁡                              (                                                      y                    -                    i                                    ,                  k                                )                                      ·                          tp2i              ⁡                              (                i                )                                                                        (        12        )            
For the outer region of the detector channels k, the following holds true:                               Low          ⁢                      -                    ⁢                      pass2            ⁡                          (                              y                ,                k                            )                                      =                              ∑                          i              =                              -                                  L                  2                                                                    L              2                                ⁢                                    Dec              ⁡                              (                                                      y                    -                    i                                    ,                  k                                )                                      ·                          tp2a              ⁡                              (                i                )                                                                        (        13        )            
where e.g. L=9
VII) Decompression
If a data compression was carried out at the beginning of the method, a data decompression now takes place by means of an inverse procedure with respect to the data compression.
VIII) Correction
The data contained in Low-pass2(y,k) represent the correction sinogram, which is now subtracted from the input sinogram. In this case, each correction projection Low-pass2(ycomp, *) is applied to n input projections.
Result[(yxe2x88x921)*N+i,k]=S[(yxe2x88x921)*N+i,k]xe2x88x92Low-pass 2 (y,k)xe2x80x83xe2x80x83(14)
where y=1,2 . . . NPRO/N and i=1,2 . . . N
The balancing method would thus be completely outlined. In order to preclude the situation where the balancing method for its part generates artifacts in addition to a satisfactory optimizing effect through definition of the thresholds in steps IV) and V), the smoothing interval in step VI) must be chosen to be largexe2x80x94usually of the order of magnitude of half a rotation.
If the length of the smoothing interval is reduced, then the method becomes critical with regard to the generation of artifacts. By way of example, if the sharp edge of an object (e.g. skull bones) runs xe2x80x9cparallelxe2x80x9d to the projection axis for a relatively long time, then this is not identified as such in step IV), as will be illustrated later with reference to the descriptive figures. Through a smoothing in accordance with step VI) in the form of averaging over a short smoothing interval, values with a high amplitude remain as a consequence along one or more adjacent channels. Despite the smoothing with a shortened smoothing interval in step VI), these values lead to considerable corrections. That becomes apparent in the image through partial rings having approximately half the length of the smoothing interval which, as it were, are drawn out or smeared out from the object structure (e.g. of the radial skull bone) in a streaky manner.
The smoothing interval of the order of magnitude of half a rotation which is defined in step VI) also has crucial disadvantages, however:
By virtue of the fact that structures which are to be assigned to objects are eliminated, in these regions of the sinogram it is no longer possible to determine the channel amplitude to be corrected. If, by way of example, on average, half of the projections are removed in the context of the differentiation in step IV), then the correction amplitude calculated in step VI) amounts to only half of the true error, since, after all, averaging was effected over 100%. A ring to be corrected thus remains after the correction with about half the amplitude and cannot, therefore, be eliminated.
It is an object of an embodiment of the present invention, therefore, to further improve the balancing method and the CT apparatus for carrying out this method.
The abovementioned object may thus be achieved by an improvement of the balancing method
Thus, in one embodiment, a method for suppressing artifacts in computer tomography raw data is proposed, which method, in different steps, determines a correction sinogram and subtracts the latter from a starting sinogram. The determination of the correction sinogram has the following steps:
high-pass filtering of a starting sinogram in the channel direction in order to filter out the long-wave structures caused by anatomical objects,
first low-pass filtering of the high-pass-filtered sinogram in the projection direction in order to improve the signal-to-noise ratio,
formation of the magnitude of a weighted gradient of each data point in the low-pass-filtered sinogram both in the projection direction and symmetrically about the corresponding channel axis, which specifies the amplitude of a change in an arbitrary direction in the sinogram, and elimination of the data point if the change amplitude thereof exceeds a first defined threshold value,
removal of residual data points in the low-pass-filtered sinogram if their amplitude exceeds a second defined threshold value,
second low-pass filtering (long-range smoothing) of the resulting sinogram in the form of averaging in the projection direction.
The correction sinogram thus obtained is finally subtracted from the starting sinogram.
This method allows an excellent identification of object structures in particular by means of the third step, so that said structures are removed in the correction sinogram and, accordingly, left in the starting sinogram, that is to say are not corrected. A further advantage is that, in the last step of the method, the smoothing interval is permitted to be significantly shortened without the method becoming critical with regard to the generation of artifacts, since, after all, all the object structures have been identified as such and removed in the correction sinogram. Reducing the smoothing interval also has the advantage, however, that the number of projections to be buffered in the computer can be drastically reduced owing to the averaging over fewerxe2x80x94to be precise significantly fewerxe2x80x94projections. The advantage consists in sparing the resources that are to be provided for the computer tomography, such as, for example, the main memory or the data transfer in the apparatus. In an advantageous manner, this also reduces the delay from the first data recording up to the generation of the CT image in the context of an xe2x80x9cimage reconstruction pipelinexe2x80x9d.
In order to obtain a better correction quality, the high-pass filtering step may advantageously be followed by median filtering of the high-pass-filtered sinogram in the channel direction.
The gradient is advantageously weighted by empirically defined scale parameters with regard to the channel or projection direction.
It is furthermore advantageous if the smoothing interval is significantly shortened in the second low-pass filtering step.
In order to reduce the computation time, as first step of the method, the number of projections used for the correction can be reduced by an averaging compression in the projection direction and can be decompressed again after the entire correction by use of an inverse procedure.
The method according to an embodiment of the invention can be implemented on a computer of a computed tomography apparatus in which the individual signal processing steps are carried out.
The method according to an embodiment of the invention can furthermore also be realized as a computer program.